function [ed, new_ed]= get_design_of_experiment

    numberOfExperiments = 20;
    numberOfParamentrs = 4;
    maxNumberOfReactions = 9;
    maxThermodinamics = 10000;
     
    reactionsCountArr = [3 5 7 9];
    transReactionsNotUsedArr = [0 1];
    thermodinamicsTreshholdArr = [5 15 50 100 10000];
    minBiomassArr = [0.2 0.4 0.6 0.8];
        
    Th = 0.0001;
    alpha = 0.14;
    imax = 2;
   % ed = get_lhs_for_each_reaction_length(numberOfExperiments / length(reactionsCountArr));
    ed = get_full_random_lhs();
    
     new_ed = sim_annealing(ed);
   
    function lhs = get_lhs_for_each_reaction_length(numberOfrows)
       lhs = [];
        for irow = 1:length(reactionsCountArr)
           lhs = [lhs; get_random_lhs_for_reaction_length(reactionsCountArr(irow), numberOfrows)];
        end    
    end    
   
   
   function lhs = get_random_lhs_for_reaction_length(reactionLength, numberOfrows)
        lhs = [];
        for irow = 1:numberOfrows
            row = [reactionLength transReactionsNotUsedArr(randi(2,1)) ...
                thermodinamicsTreshholdArr(randi(5,1))  minBiomassArr(randi(4,1))];
            
            if (isempty(lhs))
                endFlag = 0;
            else
                endFlag = 1;
            end
            
            while endFlag
                [C,ia] = setdiff(lhs,row,'rows');
                if ~isequal(length(lhs(:,1)),length(ia))
                   row = [reactionLength  transReactionsNotUsedArr(randi(2,1)) ...
                        thermodinamicsTreshholdArr(randi(5,1))  minBiomassArr(randi(4,1))];
                else
                    endFlag = 0;
                end
            end
            lhs = [lhs; row];
        end
    end 
   
    function lhs = get_full_random_lhs()
        lhs = [];
        for irow = 1:numberOfExperiments
            row = [reactionsCountArr(randi(4,1))  transReactionsNotUsedArr(randi(2,1)) ...
                thermodinamicsTreshholdArr(randi(5,1))  minBiomassArr(randi(4,1))];
            
            if (isempty(lhs))
                endFlag = 0;
            else
                endFlag = 1;
            end
            
            while endFlag
                [C,ia] = setdiff(lhs,row,'rows');
                if ~isequal(length(lhs(:,1)),length(ia))
                   row = [reactionsCountArr(randi(4,1))  transReactionsNotUsedArr(randi(2,1)) ...
                        thermodinamicsTreshholdArr(randi(5,1))  minBiomassArr(randi(4,1))];
                else
                    endFlag = 0;
                end
            end
            lhs = [lhs; row];
        end
    end
    
    function f = doe_function(reactionsNumber, transReactionsNotUsed, thermodinamicsTreshhold, minBiomass)
        f = (1 - reactionsNumber / maxNumberOfReactions) + transReactionsNotUsed + ...
            (1 - thermodinamicsTreshhold / maxThermodinamics) + minBiomass;
    end

    function Dc = sim_annealing(Dc)
        flag = 1;
        i = 1;

        while(flag == 1)
            col1 = randi(numberOfParamentrs,1);
            col2 = randi(numberOfParamentrs,1);

            while(col2 == col1)
                col2 = randi(numberOfParamentrs,1);
            end

            row1col1 = randi(numberOfExperiments,1);
            row2col1 = randi(numberOfExperiments,1);
            row1col2 = randi(numberOfExperiments,1);
            row2col2 = randi(numberOfExperiments,1);

            while(row1col1==row2col1)
                 row2col1=randi(3,1);               
            end

            while(row1col2 == row2col2)
                 row2col2 = randi(numberOfExperiments,1);    
            end

            Dnew = Dc;
            Dnew(row1col1, col1) = Dc(row2col1, col1); 
            Dnew(row2col1, col1) = Dc(row1col1, col1);
            Dnew(row1col2, col2) = Dc(row2col2, col2);
            Dnew(row2col2, col2) = Dc(row1col2, col2);

            value_h_Dc = 0;
            for j = 1 : numberOfExperiments
                value_h_Dc = value_h_Dc + doe_function(Dc(j,1), Dc(j,2),Dc(j,3),Dc(j,4));
            end 
            func_h_Dc = value_h_Dc / numberOfExperiments;
 
            value_h_Dnew = 0;
            for j = 1 : numberOfExperiments;
                value_h_Dnew = value_h_Dnew + doe_function(Dnew(j,1), Dnew(j,2),Dnew(j,3),Dnew(j,4));
            end 
            func_h_Dnew = value_h_Dnew / numberOfExperiments;
           
            
            delta_h = func_h_Dnew - func_h_Dc;
            tmprand = rand(1,1);
            tmp = exp(-delta_h / Th);

            if(delta_h <= 0)
               Dc = Dnew;
               flag = 1;
            elseif(tmp >= tmprand)
               Dnew = Dc;
               flag = 0;
            else
               i = i + 1;
            end
            if( i > imax)
               Th = alpha * Th;
               i = 0;
            end
        end
    end
end

